This document relates to constrained dynamics simulation.
In dynamics simulations, constraints are often employed on the allowable motions of the bodies being simulated. In molecular dynamics (MD) simulations, constraints may be employed to increase the size of the simulation time step without substantially impacting the accuracy of the simulation. Generally, constraints can eliminate fast vibrational degrees of freedom from a system being simulated. For example, the distance between two bonded particles is often constrained to be constant when one of the particles represents a hydrogen atom.
Ignoring the constraints, velocities and positions can be updated in a simulation using equations of motion, for example, that specify the force on each body based on its surrounding bodies. But these updated positions and velocities will not in general satisfy the constraints. One approach to introducing the constraints is to correct the positions of the bodies after an update to move the bodies to positions that would have been achieved had the constrained equations of motion been satisfied throughout the time step. An algorithm for such a correction is the SHAKE algorithm; there are other similar algorithms. In general, the SHAKE algorithm is an iterative optimization that works on a per constraint basis to move the particles incrementally to improve the matching of the constraints. Having corrected the updated positions, the computed velocities are generally then inconsistent with the positions so the velocities are updated based on the change in position in the time step. A further update of the velocities can optionally be performed to require that the velocities are essentially perpendicular to the constrained directions by removing velocity components in the constrained directions, but such a further update is not always required.
A conventional constrained NVE (constant energy) velocity Verlet dynamics simulation algorithm that takes this approach can be described by the following time step procedure (in general, in equations below the notation a_b denotes a subscript b, and the notation a^b denotes a superscript b). Let:                F(X) be a function that produces the vector of forces given the vector of positions;        X—0 and V—0 be vectors that give all the particle positions and velocities at the beginning of a timestep and xi—0 and vi—0 give the position and velocities of just particle i;        DT be the duration of the timestep;        M be the particle mass matrix (and mi the mass of particle i);        X—1 and V—1 be the positions and velocities at the end of a timestep; and        Vectors like X—1*, V_h* and so forth be computational intermediates.        
In these terms, the conventional constrained NVE velocity Verlet time step involves the following procedure:                V_h*=V—0+M^−1 F(X—0) DT/2 . . . half advance velocity        X—1*=X—0+V_h*DT . . . advance position        +X—1=CONSTRAIN POSITIONS . . . apply position constraints        +V_h=(X—1−X—0)/DT . . . correct the velocities        V—1*=V_h+M^−1 F(X—1)DT/2 . . . half advance velocity        V—1=CONSTRAIN VELOCITIES . . . apply velocity constraints        
MD simulations are often implemented with double precision arithmetic in order to reduce non-physical artifacts to an acceptable level. For example, an MD simulation of an energy conserving system usually does not exactly conserve energy due in part to finite precision arithmetic. The energy will drift as a function of simulation time and the energy drift rate may be unacceptably high if a constrained MD simulation is done with single precision arithmetic.
However, many commodity computers perform single precision arithmetic significantly faster than double precision arithmetic and have limited bandwidth between memory and processor. Significant performance benefits can be realized by using single precision arithmetic due to the faster processing and lower memory bandwidth required if the energy drift can be made acceptable on constrained MD simulations.
Many MD codes implement their constraint algorithms such that round-off errors in particle position and velocity computations become directly correlated. This in turn can cause unacceptably high negative energy drift in single precision constrained MD simulations. For example, the velocity correction step, V_h=(X—1−X—0)/DT, shown above, can cause correlation and/or high energy drift.